Microscopic derivation of transport coefficients and boundary conditions in discrete 
drift-diffusion models of weakly coupled superlattices 
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A discrete drift-diffusion model is derived from a microscopic sequential tunneling model of charge 
transport in weakly coupled superlattices provided temperatures are low or high enough. Realistic 
transport coefficients and novel contact current-field characteristic curves are calculated from micro- 
scopic expressions, knowing the design parameters of the superlattice. Boundary conditions clarify 
when possible self-sustained oscillations of the current are due to monopole or dipole recycling. 
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I. INTRODUCTION 



At present, the theory of charge transport and pattern 
formation in superlattices (SL) is in a fragmentary state. 
On the one hand, it is possible to establish a quantum 
kinetic theory from first principles by using Green func- 
tion formalisms Jl|,||. However, the resulting equations 
are hard to solve, even numerically, unless a number of 
simplifications and assumptions are made (|||. These 
include: (i) a constant electric field, (ii) simplified scat- 
tering models, and (iii) a stationary current through the 
SL. These assumptions directly exclude the description 
of electric field domains and their dynamics although im- 
portant results are still obtained ||,0J. The stationary 
current density probes the difference between strongly 
and weakly coupled SL. It also indicates when simpler 
theories yield good agreement with quantum kinetics. 
The main simpler theories are (see Figure 1 of Ref . H ) : 

(i) Semiclassical calculations of miniband transport using 
the Boltzmann transport equation || or simplifications 
thereof, such as hydrodynamic or drift-diffusion fij] 
models. These calculations hold for strongly coupled SL 
at low fields. In the miniband transport regime, elec- 
trons traverse the whole SL miniband thereby perform- 
ing Bloch oscillations and giving rise to negative differ- 
ential conductivity (NDC) for large enough electric fields 
||. The latter may cause self-sustained oscillations of 
the current due to recycling of charge dipole domains as 
in the Gunn effect of bulk n-GaAs HI). 

(ii) Wannier-Stark (WS) hopping transport in which elec- 
trons move parallel to the electric field through scatter- 
ing processes including hopping transitions between WS 
levels flUTf . Calculations in this regime hold for inter- 
mediate fields, larger than those corresponding to col- 
lisional broadening of WS levels, but lower than those 
corresponding to resonant tunneling. 

(iii) Sequential tunneling calculations valid for weakly 



coupled SL (coherence length smaller than one SL pe- 
riod) at basically any value of the electric field [p]-fi3| . 
A great advantage of this formulation as compared with 
(i), (ii) or Green function calculations is that boundary 
conditions can be derived naturally and consistently from 
microscopic models . 

On the other hand, the description of electric field do- 
mains and self-sustained oscillations in SL has been made 
by means of discrete drift models. These models use sim- 
plified forms of the tunneling current through SL barriers 
and discrete forms of the charge continuity and Poisson 
equations Although discrete drift models yield 

good descriptions of nonlinear phenomena in SL, bridging 
the gap between them and more microscopic descriptions 
p2| , p3| is greatly desirable for further advancing both the- 
ory and experiments. 

A step in this direction is attempted in the present pa- 
per. Our starting point is a microscopic description of 
a weakly coupled SL by means of discrete Poisson and 
charge continuity equations. In the latter the tunneling 
current through a barrier is a function of the electro- 
chemical potentials of adjacent wells and the potential 
drops in them and in the barrier. This function is derived 
by means of the Transfer Hamiltonian method provided 
the intersubband scattering and the tunneling time are 
much smaller than the typical dielectric relaxation time 
|]l2| . From this microscopic model and for sufficiently low 
or high temperatures, we derive discrete drift-diffusion 
(DDD) equations for the field and charge at each SL pe- 
riod and appropriate boundary conditions. The drift ve- 
locity and diffusion coefficients in the DDD equations are 
nonlinear functions of the electric field which can be cal- 
culated from first principles for any weakly coupled SL. 
These equations are of great interest for the study of non- 
linear dynamics in SL. They are simpler to study than 
microscopic model equations for which only numerical 
simulation results are available mm . 
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In the present work, natural boundary conditions for 
DDD equations are derived from microscopic calcula- 
tions for the first time. Previous authors had to pro- 
pose boundary conditions with adjustable parameters 
which gave qualitative agreement with experimental re- 
sults [ft3|-[l7ll . The present boundary conditions relate 
current density and field at contacts and can be calcu- 
lated for a given configuration of emitter and collector 
contact regions. As it is well-known, boundary condi- 
tions select the stable charge and field profiles in the 
SL, and therefore are crucial to understand which spatio- 
temporal structures will be observed in the SL for given 
values of the control parameters p~3|— |l7|| . 

The rest of the paper is as follows. In Section ||, 
we review the microscopic sequential resonant tunneling 
model. We obtain the minimal set of independent equa- 
tions and boundary conditions describing this model. 
Our derivation of the DDD model is presented in Sec- 
Numerical evaluation of velocity, diffusion and 



tion 



III 



contact coefficients for several SLs is presented in Section 
IV. Section [v] contains our conclusions. The Appendix 
contains an evaluation of the transport coefficients for 
negative values of the electric field. 



II. MICROSCOPIC SEQUENTIAL TUNNELING 
MODEL 



The main charge transport mechanism in a weakly cou- 
pled SL is sequential resonant tunneling. We shall assume 
that the macroscopic time scale of the self-sustained os- 
cillations is larger than the tunneling time (defined as 
the time an electron needs to advance from one well to 
the next one). In turn, this latter time is supposed to be 
much larger than the intersubband scattering time. This 
means that we can assume the process of tunneling across 
a barrier to be stationary, with well-defined Fermi-Dirac 
distributions at each well, which depend on the instanta- 
neous values of the electron density and potential drops. 
These densities and potentials vary only on the longer 
macroscopic time scale. 



A. Tunneling current 



The tunneling current density across each barrier in 
the SL may be approximately calculated by means of the 
Transfer Hamiltonian method. We shall only quote the 
results here |L2]]. Let eJ &i i and eJ^.c be the currents in 
the emitter and collector contacts respectively, and let 
eJ^j+i be the current through the zth barrier which sep- 
arates wells i and i + We have 
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(3) 



• i = X, . . . , N — 1, n is the number of subbands 
in each well i with energies e l c - (measured with 
respect to the common origin of potential drops: 
e = at the bottom of the emitter conduction 
band). e F = h 2 {3n 2 N D )i /(2m*) are the Fermi 
energies of the emitter and collector regions calcu- 
lated as functions of their doping density Njj. m* w 
and m* are the effective masses of the electrons at 
the wells and barriers, respectively. 



j are given by 
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where fcj and on are the wave vectors in the wells 
and the barriers, respectively, fc, depends on the 
electric potential at the center of the ith well, Wu 
whereas a,i depends on the potential at the begin- 
ning of the ith barrier, Wi + V Wi /2. See figure |j. 
Vi and V Wi , i = 1, . .. , N, are the potential drops 
at the ith barrier and well, respectively. We as- 
sume that the potential drops at barrier and wells 
are non-negative and that the electrons are singu- 
larly concentrated on a plane located at the end 
of each well (which is consistent with the choice 
of on; the choice of fcj is dictated by the Trans- 
fer Hamiltonian method). The potential drops Vq 
and Vn correspond to the barriers separating the 
SL from the emitter and collector contacts, respec- 
tively. Ai = eV Wo — 2eWo is the energy drop at 



2 



the emitter region, and eVb is the barrier height in 
the absence of potential drops. 

• Ti is the dimensionless transmision probability 
through the ith barrier separating wells i and i + 1: 



Ti(e) 



k 2 

K i+ 



(8) 



provided otid~3> 1. 



• w and d are the widths of wells and barriers respec- 
tively. 

• Scattering is included in our model by means of 
Lorentzian functions: 
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(e- £ y 2 +7 2 



(9) 



(for the ith well). The Lorentzian half- width is 
7 = h/T sc , where t sc is the lifetime associated 
to any scattering process dominant in the sample 
(interface roughness, impurity scattering, phonon 
scattering. . . ) |19|-|22||. For the samples considered 
here, 7 ranges from 1 to 10 mcV |17| . Of course this 
phcnomenological treatment of scattering could be 
improved by calculating microscopically the self- 
energy associated to one of the scattering proc- 
cesses mentioned above. However this restriction 
to one scattering mechanism would result in a loss 
of generality and simplicity of the model. 

• The integration variable e takes on values from the 
bottom of the ith well to infinity. 

Of course this model can be improved by calculat- 
ing microscopically the self-energies, which could include 
other scattering mechanisms (e.g. interface roughness, 
impurity effects j|,[l3)) or even exchange-correlation ef- 
fects (which affect the electron charge distribution in a 
self-consistent way) . We have assumed that the electrons 
at each well are in local equilibrium with Fermi energies 
e,„ , which define the electron number densities n; : 



.MrT 



A* cl (e)In 



1 



de, (10) 



Notice that the complicated dependence of the wave vec- 
tors ki and at with the potential, Wi , may be transferred 
to the Fermi energies by changing variables in the inte- 
grals of the system (|j) so that the lower limit of integra- 
tion (the bottom of the ith well) is zero: e' = e + eWi. 

Then the resulting exp ressions have the same forms as 

J+i 



Equations (f 
replaced by 



and ( p. 0D if e l ci , €q- , and e Wi in them are 



eci — e ci 



(11) 
(12) 



respectively. Wi is given by (0). The integrations now 
go from e! = to infinity. Notice that ecj is independent 
of the well index i provided we assume that the energy 
level drops half the potential drop for the whole well eV Wi 
with respect to its position in the absence of bias. Eq. 
( |l0| ) becomes 



ni(ni) 



m* w k B T 



A C i(e) In 



1 + e' fc s T 



de. (13) 



Here Aci(e) is obtained by substituting eci (the energy 
of the first subband measured from the bottom of a given 
well, therefore independent of electrostatics) instead of 
e l ci in (H). Notice that ( [l3|) defines a one-to-one relation 
between rij and ^ which is independent of the index i or 
the potential drops. The inverse function 

fM = n(m,T), 

gives the chemical potential or free energy per electron. 
This is the entropic part of the electrochemical potential 
(Fermi energy) 



i-l 



e Wi = ^{ni,T) ~ ey^Vj + V Wj ) ~ 



eV,,, 



(14) 



According to (|l4| ) , the Fermi energy, e Wi (electrochemical 
potential), is the sum of the electrostatic energy at the 
ith well, -e^jloW + Vwj) ~ e Vwi/2, and the chemical 
potential, /i; = fi(rii,T). 

After the change of variable in the integrals, the wave 
vectors in (0) become: 



tiki 



hot; 



1 2m* eV b 
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V w 
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1 2m* eV h 



eVi - eV Wi+1 - e 



(15) 



where now e = at the bottom of the ith well. This 
shows that the tunneling current density, Jj,i+i, in (j^) 
is a function of: the temperature, [ii and /Ltj+i (therefore 
of Ui and n^+i), the potential drops Vi, V^+i, V Wi , and 

Jia+i =3(ni,ni+i,Vi,Vi+i,V Wi ,V Wi+1 ). (16) 
Similarly, we have 

J M =S e (m,JVb,V r ,V r Iol ) ) (17) 
JN, c = ^c{n N ,N D ,V N ,V WN ). (18) 
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B. Balance and Poisson equations 

The 2D electron densities evolve according to the fol- 
lowing rate equations: 



~~dt 



Ji- 



J, 



1. 



,N. (19) 



The voltage drops through the structure are calculated as 
follows. The Poisson equation yields the potential drops 
in the barriers, Vi, and the wells, V W i (see Fig. pi): 



Vi-i , e{n t -N™) 



w £ d ' 2 
V = Vj-i | e(m-N%) 
d d e 



(20) 
(21) 



where e w and e are the GaAs and AlAs static permit- 
tivities respectively, rii is the 2D (areal) electron number 
density (to be determined) which is singularly concen- 
trated on a plane located at the end of the ith well, and 
Nf, is the 2D intentional doping at the wells. 



C. Boundary conditions 

The emitter and collector layers can be described by 
the following equations: 

e^Ai seV V 

— — = —j- , cr e = 2e — ~ eN(e F )Aidi , (22) 



e w A 2 eV N eN D S 2 s w e F 



eS 2 



eS 3 



e F ( 1 

a c = 2e w — = eN D S 2 + -63 
eo3 V 2 



(23) 
(24) 



To write the emitter equations ( p2| ) , we assume that there 
are no charges in the emitter barrier B3| . Then the elec- 
tric field across Si (see Fig. [j]) is equal to that in the emit- 
ter barrier. Furthermore, the areal charge density a e re- 
quired to create this electric field is provided by the emit- 
ter. N(e F ) = m^U^ 2 (3N p> / -k a )z is the density of states 
at the emitter Fermi energy e F = % 2 {?>it 2 Nd) 3 / {2m*,). 
The collector equations (|2^) and ( f2~4| ) ensure that the 
electrons tunneling through the iVth (last) barrier are 
captured by the collector. They hold if the bias is large 
enough (see below). We assume that: (i) the region of 
length 8 2 in the collector is completely depleted of elec- 
trons, (ii) there is local charge neutrality in the region of 
length 8 3 between the end of the depletion layer S 2 and 
the collector, and (iii) the areal charge density a c required 
to create the local electric field is supplied by the collec- 
tor. Notice that eN D (S 2 + |<5 3 ) in is the positive 2D 
charge density depleted in the collector region. Equa- 
tions (|23|) and (p4[) hold provided Vn > e w e F d/ '(eeS 3 ), 



A 2 > 0, S 2 > and ^3 > 0. For smaller biases resulting 
in Vn < s w e F d/(eeS3), a boundary condition similar to 
(p2|) should be used instead of (23) and (§3): 



e w A 2 eeV N V N , >r 5 , . 

— = — = — - 1— , <J C = 2e — ~ eN(e F )A 2 5 2 . (25) 
02 d 



Notice that A 2 and S 2 have different meanings from A 2 
and S 2 in (p^). 

The condition of overall voltage bias closes the set of 
equations: 



N N 



v = ^v i + ^v w 

i=Q i=l 



Ai + A 2 + e F 



(26) 



This condition holds only if Vn > £ w c.Fd/(ee5$)\ other- 
wise (A 2 + e F ) should be replaced by A 2 in (]26]). 

Notice that we can find 8% and Ai as functions of Vo 
from (||): 

Ai = = Vq, Si indetermined or 



5i = J^M^= hn T £w) \ , Al = ^.(27) 



e 2 N(e F ) em *J(3 NlD )i 



Enid 



Similarly we can find S3 by solving (|23|) and ( |24| ) in terms 
of Vjv arid Np . From this equation and ( p3[ ) , we can find 
S 2 and A 2 as functions of Vn- 



2s 



eN D d 



2 2e w e F N D d 2 



1} 



<5 2 = 



A 2 = 



e 2 N D S 3 \ s w e F d 
2ee 2 F ( e£V N S 3 
e 2 Nr)8 2 \ e w e F d 



- 1 



- 1 



/ eeV N S 3 
V £ w £Fd 
( eeV N S 3 
\ e w e F d 



N 



- 1 



(28) 



where 6(x) is the Heaviside unit step function. The 

boundary conditions ( |23| ) and (|24|) do not hold 

if eeVN&s < e w e F d. This occurs if Vn < 
1 5 

(37r 2 )3^M7e7;/[2eV2^<]- In this case, we should 
impose the alternative boundary conditions (p5|). From 
these, we obtain 



S 2 = 



2e„ 



Hits (2s^ 



e 2 N(e F ) 



A 2 = 



eeVN$ 2 



em-i 



1 - 



e w e F d 



(29) 



The critical potential Vjv = £ w £ Fd/ (ee5 3 ) corresponds to 
Vn = £u)£f'^/(v / 3 ee5 2 ). There is a small mismatch be- 
tween ( p^ ) and ( p9| ) at this critical potential: eeVjv /d = 
e w e F /S3, but eeVN /d ^ s w t F jS 2 . This imperfection can 
be fixed by using a more precise relation between the 
charge at the collector cr c , A 2 and S 2 , but we choose not 
to delve more in these details. In all cases, we have shown 
that the potential drops at the barriers separating the SL 
from the contact regions uniquely determine the contact 
electrostatics. 
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In Ref. O] global charge conservation: 

N 1 
<r e + 5^(ni - eA^) = eJVr,(<J 2 + -S a ) , (30) 

»=i 

was used instead of (^J) [which is a condition similar to 
the one we impose at the emitter contact, fl22])]. Sub- 
stitution of (^J) instead of ( |30| ) modifies minimally the 
numerical results reported in Refs. |L2| and |l8| ]. 

D. Elimination of the potential drops at the wells 



Notice that the three last equations are constitutive re- 
lations obtained by substituting (^lj) in the functions S, 
S e and S c of (|l6|), ( |l7j ) and (|l8|), respectively. The func- 
tions Ai(Vb) and A 2 (V/v) are given by ( pjj ) and (p8|), re- 
spectively. Equations (B3h for i = 0, N may be considered 
the real boundary conditions for the barriers separating 
the SL from the contacts. These boundary conditions are 
the balance of current density including special tunneling 
current constitutive relations J ei i and Jn,c- The latter 
depend on the electron densities at the extreme wells of 
the SL and the potential drops at the adjacent barriers. 



The previous model has too many equations. We can 
eliminate the potential drops at the wells from the sys- 
tem. For (|20|) and (§l|) imply 

e w V Wi Vi-i + V, 



ew 

Then the bias condition (]2( 



V= 1 



ew 
e w d 



N 

E 

i=0 



V, 



2d 

) becomes 

£ {V + V N ) w 



(31) 



2e w d 
A 2 + 



(32) 



where Ai = Ai(Vq) and A 2 = A 2 (V/v)- Instead of the 
rate equations (|19|), we can derive a form of Ampere's law 
which explicitly contains the total current density J(t). 
We differentiate (|l]) with respect to time and eliminate 
rii by using ([Ts)) . The result is 



e dV, 
ed dt 



J, 



i = Q,l,...,N, (33) 



where eJ(t) is the sum of displacement and tunneling cur- 
rents. The time-dependent model consists of the 2N + 2 
equations (^l|) , (32) and (|3^ ) [the currents are given by 
Eqs. (|), ([H]), @, (H) and ©], which contain the 
2N + 2 unknowns m (j = 1, . . . , N), Vj (j = 0, 1, . . . , N), 
and J. Thus we have a system of equations which, 
together with appropriate initial conditions, determine 
completely and self-consistently our problem. For con- 
venience, let us list again the minimal set of equations 
we need to solve in order to determine completely all the 
unknowns: 



e dVi 
ed dt 

W 
d 



j, 



J(t), 



i-1 



e (fa — N]j) 



e 



V= 1 



ew 
+ 



N 

E 

;=o 



Vi- 



-0,1, 

i = l, 
(Vb 



Vn) ew 



Ai(^o) 



2e w d 
A 2 (V N ) + 



e 

, n,4-i , Vs- 



H,i+1 = ^K'Hi 'H+l, fi-li Vi, Vj+i), 

J e ,i = 3 e (ni, Vb,Vi), 
Jn,c — 2 c (njv, Vjv_i, Viv). 



(34) 



(35) 

(36) 
(37) 
(38) 



III. DERIVATION OF THE DISCRETE 
DRIFT-DIFFUSION MODEL 

It is interesting to consider the relation ( |l0| ) between 
the chemical potential and the electron density at a well 
for different temperature ranges: 



fa(fa 



m* w k B T 



A C i(e) In 



1 + e 1 



de. 



Assuming that fa 3> k B T, we may approximate this ex- 
pression by 



ni(pi) 



7r 2 h 2 Jo 



Aci{e) (fa - e) de. 



(39) 



Thus fa approaches a linear function of fa if fa 3> k B T. 
For the SL used in the experiments we have been refer- 
ring to, fa — e is typically about 20 meV or 232 K. Thus 
"low temperature" can be "high enough temperature" in 
practice. Provided the Lorentzian Aci(e) is sufficiently 
narrow, Aci(e) ~ irS(e — ecx), so that 



nh 2 7 



if (fa - eo) > k B T, eci > 7- (40) 



Interestingly enough, a linear relation between fa and fa 
also holds at high temperatures. To derive it, notice that 
l n(l + e x ) ~ In 2 + x/2 if x <C 1 and use this relation in 
&■ 

n i(p-i) ~ 7r~5T2 / A C i(e)(2k B Tln2 + fa-e)de. (41) 

271~ 2 ft Jo 



If we now set A C \(e) ~ ir S(e — eci), the result is 
-2fc s Tln2 



2'Kfi 2 rij 



fa 



if (fa - eci) < k B T < (e C2 - eci), and e C i > 7- Tne 
additional condition (thermal energy small compared to 
the difference between the energies of the two lowest sub- 
bands) is needed to keep all electrons in the first subband. 
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For otherwise the second subband may be populated and 
Equation (|l3| ) should be transformed accordingly. Thus 
our "high temperature" approximation can be satisfied 
in SL with large enough energy differences (ec2 — eci)- 

A different approximation is obtained if we first impose 
that eci > 7: 



m* w k B T 
riiifii) w — — In 



This yields 



Mi 



eci + k B T In 



and therefore 

Mi - eci ? 



Trh 2 n; 



m* w kBT 



if % m^> m*k B T, 



if 7i rij <C m* w k B T. 



At low temperatures, the chemical potential again de- 
pends linearly on the electron density according to (40), 
whereas it has ideal gas logarithmic dependence at high 
temperatures. 

The same considerations used to obtain (^) or ( Ell ) 
would indicate that the electron flux across the ith bar- 
rier becomes 



n 2 v. 



(/) 



n i+ iv. 



W 



,(/) 



either at low or high enough temperatures. Here v\ 

and are functions of Vi, Vi±±. They have dimen- 
sions of velocity and correspond to the forward and back- 
ward tunneling currents which were invoked in the deriva- 
tion of phenomenological discrete drift models. When 
£», = e Wi+1 , or equivalently, /Ltj+i = Mi + eV i + e (Ku + 
V Wi+1 )/2, — according to (g). Equation jl^ ) 

implies that Mi+i = Mi if n i+i = n i, and therefore 
we conclude that v^' — at zero potential drops 
Vi + (V Wi + V Wi+1 )/2 = 0. Notice that e Wi+1 — e becomes 
Mi+i ~ eVi — eV Wi+1 /2 — e' after changing variables in 

the integral (||). Then v+ is approximately zero unless 
< e C i + eV w j2 < fi i+1 - eV l - eV Wi+1 /2. For volt- 
ages larger than those in the first plateau of the current- 
voltage characteristic curve this condition does not hold. 
In fact for these voltages, the level CI of well i is at a 
higher or equal potential than the level C2 of well i + 
Then e C i > Mi+i ~ eVi ~ e(V Wi + V Wi+1 )/2. 

The previous results yield DDD models with the po- 
tential drops at the barriers and the total current density 
as unknowns, the same as in Eqs. (^) - (^8|). The main 
difference with previously used discrete drift models is 
that the velocity depends on more than one potential 
drop. To obtain these simpler models, we further assume 
that eVi/ed and eVi±\/ed are approximately equal to an 



average field Fj (e is an average permittivity to be cho- 
sen later). Then V Wi — weFi/e w according to (j3l|). This 
assumption departs from previous approximations and 
yields a new model. The point of contact with our previ- 
ous results is that Aci{e)Acj(e + eVi + e\V Wi + V Wi+1 }/2) 
is the controlling factor in the expressions for c") and 
i;W (the transmission coefficient contains an exponential 
factor, e _2Qid , which is almost constant at the energies 
contributing most to the integral). This controlling fac- 
tor is uniquely determined by the potential drop 



Vi 



v w 



V 



W 

c w 



sFi = {w + d)Fi, 



provided we define the average permittivity as 
d + w 



d 1 w 



(42) 



This expression corresponds to the equivalent capaci- 
tance of two capacitors in series. Thus the behavior of 
forward and backward drift velocities is most influenced 
by the potential drop Vi + (V m + V Wi+1 )/2 w Fj(d + w) 
and the new DDD model (see below) should yield results 
similar to those of the microscopic sequential tunneling 
model. We have 



M+i 



d + w 
_ riiv(Fi) n i+ i - n t 



- , , , ,2 D(Fi), (43) 

v(F) = v if) (F) -v^(F), D(F) = (d + w)v ib \F). (44) 

To calculate v^(F) and v^(F) from J i<i+1 in (g), we re- 
place e Wi , e* ci , e Wi+1 and e^ 1 by Mi, eci, [ii+i-e(d+w)F 
and ecj — e(d + w)F, respectively. The wavevectors in 
the integrand should be 

Tiki = y / 2m* u e 



Hon 



2m [eV b -^—-e 



hki+i = ^2m* u [e + e(d + w)F] 



hcti-i 



2m* 



eV b + e(d+^)F-e 



frcii+1 — \ 2m* 



(45) 



and the integration variable e ranges from to 00. We 
substitute Mi( n i) according to ( |l3| ) in the result. Then 
we obtain a function J{rii, nj+i, F): 



J(ni,n i+1 ,F) = 5 ni,m+i, 



eFd eFd eFd\ 



e ) 



(46) 



G 



[equivalent to setting Vi = eFd/e, or Vi + (V Wi + 
V Wi+1 )/2 — (d + w)F after transforming this formula 
to the form (pq)]. Notice that (as said above) 

„</)(()) =^(0) = ^- 



for the sake of simplicity. For fields higher than the first 
resonance, D(F) « 0, and ( |50|) becomes 



s dF, 
e dt 



niv(Fi] 



= m , 



(53) 



for the tunneling current to vanish at zero field and equal 
electron densities at adjacent wells. Furthermore, no- 
tice that D(F) vanishes if eci — Hi+i > —e + V Wf + 
V Wt+1 )/2 « -e(d + w)Fi. Thus according to @, D(F) 
vanishes if fi. 2 ni + i < m* e(d + w)Fi, which is certainly 
satisfied for all average fields larger than the first reso- 
nant field (ec2 — eci)/[e(d + w)]. In the low temperature 
limit (or in the high temperature limit mentioned earlier 
in this Section, provided it exists), we have 



J(ni,n i+ i,F) = 



■v(F) 



7i.,: + 1 - n 



(d + w) : 



-D(F). (47) 



Then we may use 
v(F) = 



D(F) 



(d + w)J(N%,N%,F) 



N 



D 



(d + w) 2 J(0,N%,F) 
NVt 



(48) 
(49) 



to calculate the drift velocity and the diffusion coeffi- 
cient from the tunneling current. The integrals from (j^) 
appearing in these expressions may be approximated by 
means of the Laplace method: we should just expand 
their controlling factor mentioned before about its maxi- 
mum value e = e(F). The resulting formulas are cumber- 
some and we choose not to write them here. We show in 
the Appendix that v^(-F) = v^(F) = D(F)/(d + w), 
and v(-F) = -v(F). 

Equations ( ffTj ) to ( ff9|) may be used in ([33]) to write 
the Ampere law as 



e dF\ njvjFi) n i+ i - n, 

e dt d + w { l) (d + w) 2 



J(t). 



(50) 



for i = 1, . . . , N — 1. Poisson equation ( J34| ) becomes 
Fi-Fi_ x = -{m-N%), 

£ 



(51) 



for i = 1, . . . , N. Equations ( |50|) and ( pl[) constitute a 
DDD model which may be analyzed on its own together 
with appropriate bias and boundary conditions. As bias 
condition we adopt 



JV 

{d + w) Y, F i = v - 

i=l 



(52) 



Notice that potential drops at the contacts are repre- 
sented only by the term Fjv (d + w). Equation (]52) is 
obtained by inserting Vi + (V Wi + V Wi+1 )/2 = (w + d)Fi 
into (p6|), and omitting 



(d + w)F a + 



Ai + A 2 + 2e F 

2e 



which is the usual discrete drift model used in previous 
theoretical studies (l^-[l7]]. 

In Section 2.1 of Ref. [pi, A. Wacker derived a for- 
mula similar to ( ^3| ) with = and v^(F) cx 
T/[(eF(d + w) + eci - ^Cjf + T 2 ], for fields compara- 
ble to (ecj — eci)/[e(d + w)]. At low fields, the res- 
onant tunneling current between levels CI of adjacent 
fields was found to be proportional to W(F) = eF(d + 
w)/[(e 2 F 2 (d+w) 2 +T 2 ] and independent of rij. While the 
first approximation of Wacker's (for fields close to higher 
resonances, CI — » Cj, j = 2, 3, . . .) is compatible with 



our result (43), the second approximation is an artifact 
of the extra unnecessary assumption e Wi — e Wi+1 p3[. 
We shall show in Section [lV| that our drift velocity ( Eg ) 
may have at low fields the same shape as function W(F) 
for certain SL; see Fig. ||(a). Then the corresponding 
stationary current for a space homogeneous field profile 
with rii — Np (which implies equality of chemical po- 
tentials at adjacent fields) will be proportional to W(F) 
and our result will agree with Wacker's (for this special 
case). Fig. |2|(b) shows that things may be different for a 
different SL configuration. 

The boundary conditions for Fq and Fjy are also 
Ampere's law but now the tunneling currents (Q) and 
(H) (from the emitter and to the collector, respectively) 
have to be used instead of (^|) . The same approximations 
as before yield 



J e ,i = ^(nx^Fod/E^Fod/E) 



;(/) 



d + w 

Jat jC = Z c (n N ,£F N d/£,£~F N d/£) 



d + w 



.,(/) 



(F N ). 



(54) 
(55) 



Notice that there is no backward tunneling from the col- 
lector region to the SL because we are assuming that the 
potential drop Vjv is larger than £ w epd/ (eec^). Assuming 
now that (pi) and (B5|) are identities, we find 



jP(F) = E, (o 
w^(F) = 



EFd EFd\ 



N 



D 



& f \F)-~AN 



EFd 



EFd\ 



N w cy D , _ 



EFd EFd\ 



(56) 
, (57) 
(58) 



The tunneling current across a barrier is zero if the Fermi 
energies of the adjacent wells are equal. This occurs if 
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the electron density at the first well takes on an appropri- 
ate value n™ such that the corresponding Fermi energy 
equals that of the emitter. Then 



and therefore 



H e «,0,0)=0, 



m nj>W(0) 
Je { 1 d + w 



IV. NUMERICAL CALCULATION OF DRIFT 
VELOCITY AND DIFFUSION 

In this Section, we shall calculate the functions v(F), 

D(F), ji f) (F), w {b) (F) and w^(F) for different SL 
used in experiments p7| . Fig. |2|(a) depicts the elec- 
tron drift velocity v(FJ for the 9nm/4nm GaAs/AlAs 
SL (9/4 SL) of Rcf. @ calculated by means of © (at 
zero temperature; m* = m* w for simplilcity) . The inset 
compares v(F) to the backward and forward velocities 
given by v^ b \F) = D(F)/(d + w) [D(F) as in @] and 
y(/)(F) = + u( fc )(F). The rapidly decreasing dif- 

fusivity D(F) determines the position and height of the 
first peak in v(F). Notice that v(F) behaves as expected 
from general considerations: it increases linearly for low 
electric fields, it reaches a maximum and then decays be- 
fore the influence of the second resonance is felt. If D(F) 
decays faster, a rather different v(F) is found. Fig. ||(b) 
shows v(F) for the 13.3/2.7 SL: there is a wide region 
before the first peak in which v"(F) > 0. 

Figures || and || show the boundary functions ji (F) , 
w^(F) and wWjF) for the 9/4 and 13.3/2.7 SL, respec- 
tively. Again they behave as expected: (i) the emitter 
forward current peaks at the resonant values of the elec- 
tric field [different from those of v^'(F)], (ii) the emit- 
ter backward tunnel velocity decreases rapidly with field, 
and (iii) the collector forward velocity increases mono- 
tonically with field given the large difference between the 
Fermi energies of the last well and the collector. 

The emitter forward current, je (F), is compared in 
Figs. | and | to the drift current, N%v(F)/(d + w), for 
different emitter doping values. Notice that the emitter 
current is systematically higher than the drift current for 
large emitter doping at positive electric fields. However, 
the total current density should remain between the first 
maximum and the minimum of the drift current. This 
means that the contact field Fb should be negative, so 
that j { e f) (F Q ) - m w^(F )/(d + to) « J, with m > N%. 
Because of Poisson equation, (|34|), F\ is larger than Fq 
and, typically becomes positive. The electric field in 
the SL increases with distance from the emitter and a 
charge accumulation layer is formed (see Figure 5 of Ref. 
[ p2[ for a similar stationary field profile corresponding 
to the full microscopic sequential tunneling model) . Self- 
consistent current oscillations in this situation will be due 
to monopole recycling Ea] . Notice that previous work on 



discrete drift models assumed a fixed excess of electrons 
in the first SL well as boundary condition [|l6],[l7) . Again 
an emitter accumulation layer appeared and monopole 
self-oscillation resulted. 

( f ) 

For smaller emitter doping, j e (F) intersects 
Npv(F) I \d + w) on its second branch, and a charge de- 
pletion layer may be formed in the SL. Then there may 
be self-oscillations due to dipole recycling. These findings 
are fully consistent with the numerical results reported 
in Ref. fj| for the 13.3/2.7 SL. That paper reported 
coexistence and bistability of monopole and dipole self- 
oscillations for the first time. Coexistence and bistabil- 
ity were found for an intermediate emitter doping range 
(crossover range) lower than those used in experiments 
A different way to obtain dipole self-oscillations is 
to decrease the well width without changing contact dop- 
ing. In this way, we have numerically checked that dipole 
self-oscillations are possible with emitter doping similar 
to those used in current experimental setups |17j . 

For the usual drift-diffusion model of the Gunn effect 
in bulk n-GaAs, the effect of boundary conditions on 
the self-oscillations of the current has been well-studied 
p4| , p5| . In particular, asymptotic and numerical calcula- 
tions for realistic metal-semiconductor contacts were per- 
formed some time ago pEfl . Despite the different equa- 
tions used in bulk semiconductors, these calculations pro- 
vide results consistent with our present findings in SL: 
a boundary condition which yields accumulation (deple- 
tion) layer near injecting contact may yield current self- 
oscillations due to monopole (dipole) recycling p5| , p4[ . 
However these similarities between discrete (SL) and con- 
tinuous (bulk) drift-diffusion models should not tempt us 
into reaching hasty conclusions: discrete and continuous 
drift-diffusion models may have rather different traveling 
wave solutions pq| . In fact, it has been shown that (de- 
pending on current and doping), the DDD model may 
have monopole wave solutions which travel in the same 
direction as the motion of electrons, in the opposite direc- 
tion, or remain stationary. In the continuum limit (con- 
tinuous drift-diffusion model), wavefronts travel always 
in the same direction as the electrons j2(| . These features 
of the DDD equations may have experimentally observ- 
able consequences which will be explored elsewhere. 



V. CONCLUSIONS 



Starting from a microscopic sequential tunneling model 
of transport in weakly coupled SL, a DDD model is de- 
rived in the limits of low or high temperature. Realis- 
tic transport coefficients and novel contact current -field 
characteristic curves are calculated from microscopic ex- 
pressions, knowing the design parameters of the superlat- 
tice. Boundary conditions select stable spatio-temporal 
charge or field profiles in the SL. In prticular, they clar- 
ify when possible self-sustained oscillations of the current 
are due to monopole or dipole recycling. 
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APPENDIX A: MODELS FOR NEGATIVE BIAS 

When a negative voltage is applied, we should make 
sure that our formulas transform appropriately. For neg- 
ative bias, the charge will be singularly concentrated on 
planes located at the beginning of the wells. Then we 
should write 



\ 



2m* 



iV b -eJ2{Vj+V Wj )-e 

3=0 



instead of (g) in the expressions (g) ■ The change of vari- 
able e' = e + eWi+i (i.e., e' = corresponds to zero 
energy at the bottom of well i + 1) in the integral (||), 
then changes the wavevectors to 

hki+i = ^2m* w e, 



hen = \ / 2m* I eVb + 



eV, 



Wi+l 



hki = \ / 2m* \ e - eVi- e 



ThOi-i = \ / 2m* eV b + eV Wi + eV l + 



Wi+l 



frui+i = \ 2m* eVb - eV i+ i 



(Al) 



instead of (p"5|). 

Given the new location of the singular charge planes 
(at the beginning of wells), © still holds, but (§l]) 
should be replaced by 



W W £,„ 



Then we find 



eVj = V UH + Vw i+1 
£,„d 2w 



(A2) 



(A3) 



instead of (|31j). Inserting this equation in the functions 
S (tunneling current under negative bias), we obtain 



new functions S*(n,,nj+i, V Wi , V Wi+1 , V Wi+2 ), instead of 
S(rii, ^i+ii Vi-i, Vi, Vi+i) valid for positive voltage. To 
obtain a reduced model, we now set 

£ci = £ ci + e(d + w)F , 



,t+l 



/i t + e(d + w)F . 



fik i+ i = v / 2m* u e, 



hat 



2m* eV b + 



ewF 



yj2m* w [e-e(d + w)F], 



hati-i — \ 2m* 



Ha i+ i =i/2m 



eV b 



F-e 



(A4) 



in the integrals (0) and let the variable of integra- 
tion e range from to oo. This is equivalent to 
setting V Wi , V Wi+1 and V Wi+2 equal to ewF/e w in 
E*(ni,n i+1 ,V Wi ,V Wi+1 ,V Wi+2 ). Equations (|), ©, @) 
and the previous definitions in this Appendix imply 



eFd eFd eFd\ 

a I TH,ni+i, , , 

see) 



eFw eFw eFw 



a n i+ i,ni, ■ 



(A5) 



The Poisson equation (A2) still yields (51). Notice that 
the symmetry (A.5) implies 



D(F) 



v(-F) = —v(F). (A6) 



Given the difference between the states at the contact 
regions and the wells, the previous arguments cannot be 
used to extend the contact current-field characteristic 
curves to negative fields. Direct calculation of (|5^) - ( |5^ ) 
by means of (1) and (||) yields the results depicted in 
Figures |^ and 4 . 
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FIG. 1. Sketch of the electrostatic potential profile in a SL. 
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FIG. 2. (a) Electron drift velocity v(F) for the 9/4 SL. In- 
set: comparison of the drift velocity (continuous line) with 
the forward (dashed line) and backward (dot-dashed line) ve- 
locities, (b) The same for the 13.3/2.7 SL. Notice that the 
backward velocity or, equivalently the diffusivity, decreases 
with electric field much more rapidly for this SL. 
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FIG. 3. Functions of the electric field appearing in the 
boundary conditions for the 9/4 SL with a contact doping 
N D = 2 x 10 18 cm -3 , (a) eji f) (F) and (b) u> (b) (F) for the 
emitter and (c) w^'(F) for the collector. 
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FIG. 5. Comparison of the drift tunneling current density, 
eNp v(F)/(d + w) (continous lines) with the emitter current 
density eji (F) (dashed lines) for the 9/4 SL with two differ- 
ent emitter dopings: (a) Nn — 2 x 10 18 cm -3 corresponding 
to monopole recycling, and (b) Nr> = 2 x 10 17 cm~ 3 corre- 
sponding to dipole recycling. 
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FIG. 4. Same functions as in Figure J for the 13.3/2.7 SL 
with a contact doping Nr> — 2 x 10 cm -3 . Notice that 
eji (F) is an increasing function since ef > (ec2 — eci) in 
this SL. 
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FIG. 6. Same functions as in Figure ^ for the 13.3/2.7 
SL (a) No — 2 x 10 18 cm -3 (monopole recycling), and (b) 
No = 10 16 cm -3 (dipole recycling). 
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